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Abstract 

A cohesive zone model (CZM) for material fatigue analysis 
is presented. This cyclic CZM contains a damage variable 
and a corresponding damage evolution law. Considering the 
influence of the constraint condition for CZM, the triaxiality 
factor is also involved in this model. The new cyclic CZM 
can describe the different damage evolution mechanisms in 
the fatigue crack initiation and fatigue crack growth process. 
A unique set of model parameters is used. The material 
damage parameters included in the damage evolution law 
can be obtained from the material strain-life curve. Then 
they are used for fatigue crack growth calculation. The 
simulation results are compared with the experimental data 
in the very low cycle fatigue regime. 
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Introduction 

In reality, numerous engineering structures and 
components are operated under cyclic loading. Fatigue 
fracture of materials and structures is a failure 
phenomenon which can cause catastrophes. Therefore, 
the fatigue life prediction and fatigue crack growth 
analysis are very important issues for material science 
and engineering. Fracture mechanics concepts can be 
used to solve such problems but modeling expenses 
might become high. Alternative approaches need to be 
introduced in this field, and the selected method 
should be simple enough in theoretical understanding 
and practical application. Synthesizing all the 
consideration, the application of the cohesive zone 
models in fatigue calculation is a promising 
development. 

Cohesive zone models (CZM) describe the material 
damage in a narrow strip zone. Originally, the idea of 
the CZM was proposed by Dugdale (1960) and 
Barenblatt (1962). In 1976, the CZM was firstly applied 
by Hillerborg (1976) in finite element calculations for 
concrete material. From then on, the CZM was used 



widely in the finite element simulation for material 
fracture. It was applied by Needleman (1987, 1990) for 
modeling inclusion debonding and interface 
decohesion in microscopic and macroscopic scales, 
respectively. Subsequently, it was employed by 
Tvergaard and Hutchinson (1992) for crack growth 
calculation in elastic-plastic solids. In the last two 
decades, the CZM was improved by many researchers 
and scientific groups (Chaboche et al. 1997, Yuan et al. 
1996, Cornec et al. 2003, Scheider et al. 2003, Scheider 
et al. 2006, Mohammed et al. 2000, Chandra et al. 2002, 
Chen et al. 2003, Song et al. 2006). Now it is a robust 
tool for material fracture analysis. 

However, a remarkable fact is that all these CZMs 
were just used for material fracture under monotonic 
loading. In the late 1990's, the CZM was firstly 
extended by De-Andres et al. (1999) for fatigue crack 
growth simulation. A damage variable was introduced 
in this CZM to describe the cohesive element 
degradation. But the formulation did not contain a 
description about the damage evolution. Soon 
afterwards, Yang et al. (2001) developed a CZM based 
on a line spring model to simulate fatigue crack 
growth. The damage variable and the corresponding 
damage evolution were intrinsically involved in the 
fatigue process. However, in this model the unloading 
and reloading stiffness had different definitions. These 
models were not convenient and accurate for 
application. Then the cyclic CZM proposed by Roe 
and Siegmund (2003) explicitly introduced a damage 
variable and a damage evolution equation. The 
stiffness and maximum cohesive traction degraded 
with the damage variable. The unloading and 
reloading path used the same definition. The 
numerical calculation results of this model have 
reproduced many basic characteristics which are 
similar to the experimental phenomena in fatigue 
fracture. Nevertheless, there was no experimental 
validation for the damage evolution process and 
various other important factors for CZM were not 
considered, such as the influences of surrounding 



31 



www.seipub.org/sas 



Solids and Structures (SAS) Volume 2 Issue 3, September 2013 



continuum element, element friction and constraint 
condition. Based on Roe and Siegmund's (2003) 
method, some authors published further 
improvements and extensions (Abdul-Baqi et al. 2005, 
Bouvard et al. 2009, Xu et al. 2009). In addition, others 
also proposed various damage evolution hypotheses 
of CZMs for fatigue simulation (Ural et al. 2009, 
Pirondi et al. 2011, Turon et al. 2007, Munoz et al. 2006, 
Beaurepaire et al. 2011). Recently, Jha and Banerjee 
(2012) presented a CZM for fatigue failure in complex 
stress states. The description of the damage and 
damage evolution were still from Roe and Siegmund's 
method (2003). This new CZM can consider the 
influence of various triaxiality conditions. But in this 
model, triaxiality was not involved in the damage 
evolution law. The work just focused on fatigue life 
analysis. Fatigue crack growth has not been calculated. 
Furthermore, the simulation results just gave some 
fatigue fracture trends, but no comparison with 
experimental results was included. 

Cohesive Zone Model Approach 

Fundamentals 

The CZM intends to describe the real physical fracture 
process by phenomenological equations. The CZM can 
be adopted for a material without or with a crack. The 
total material damage is supposed to occur in the 
cohesive zone, the continuum elements in the vicinity 
are damage free. In practice, the CZM can be applied 
for any material, any fracture mode or any model 
dimension. In the CZM, there is a constitutive 
relationship between cohesive traction T and material 
separation 5, called traction-separation law (TSL). The 
TSL represents the material deterioration occurring in 
the interface under the monotonic loading process. In 
practical application, several TSLs can be chosen, for 
example, the cubic shape (Needleman 1987), the 
exponential shape (Needleman 1990), the trilinear 
shape (Tvergaard and Hutchinson 1992) and the partly 
constant shape (Scheider et al. 2003). All the TSLs 
contain two parameters, the maximum traction To and 
the critical separation <5o. The area under the TSL 
represents the cohesive energy Fo. It can be derived as 

r n = [ d °T(s)ds 

Jo v r . If the shape of the TSL is selected, the 
cohesive parameters should be determined from 
experiments. Hence the cohesive parameters are both 
model and material parameters. In this paper, the TSL 
from Scheider et al. (2003) is taken as a basis (Fig. 1). 
More details about applying this CZM can be found in 
references (Scheider et al. 2003, Scheider et al. 2006). 




FIG. 1 TRACTION SEPARATION LAW ACCORDING TO 
SCHEIDER ET AL. (2003), GIVEN BY EQUATION (1) 

Advanced Property - Triaxiality Dependent Behaviour 

Besides using the CZM in a phenomenological way, it 
is possible to find a relationship between the 
micromechanics of a fracture process and the TSL. In 
micromechanics, the material ductile fracture is 
represented by void nucleation, growth and final 
coalescence. Based on this process, a micromechanical 
model for material ductile fracture was proposed by 
Gurson (1977), later extended by Tvergaard and 
Needleman (1984), which is now called GTN model. 
Using the GTN model to simulate a tension test for 
one cell element, a stress displacement curve can be 
derived. Such a curve can be directly used as a TSL, or 
can be fitted to another common shape, which might 
be implemented in a finite element code. For this 
micromechanical traction separation curve, some 
investigators found that it was strongly influenced by 
the triaxiality conditions. Here, the triaxiality value " 

means the ratio between the hydrostatic stress * and 

the von Mises equivalent stress °' Mises . 

°Mises (2) 

Siegmund and Brocks were the first researchers to 
propose such a triaxiality dependent CZM (1999, 2000). 
More investigations on triaxiality dependent 
behaviour were published by Anvari et al. (2006), 
Scheider (2009) and Banerjee et al. (2009). The research 
conclusions of all authors can be summarized: The 
cohesive parameters To, <5o and To are not pure material 
constants but dependent on the triaxiality conditions. 
With increasing triaxiality, the maximum traction To 
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increases, but the cohesive energy To (or critical 
separation So) decreases. 

When the triaxiality dependent CZM is applied in 
finite element simulations, another important issue is 
how to get the triaxiality value for a cohesive element. 
In cohesive elements, there are no definitions about 
hydrostatic stress and von Mises equivalent stress. So 
the triaxiality can only be obtained from the adjacent 
continuum element and has to be transferred to the 
cohesive element. 

For describing the triaxiality dependent behavior, 
equation (3) is used. This equation is taken from the 
literature (Scheider 2009) and a small modification is 
introduced. 



TJa Y 



:3.3-4.8exp 

_ 1 
~ 5 



0.1 -H 
0.45 

0.138 + 0.635 exp 



- H 



f 0.858 
[ 0.788 



+ 0.5346 



23 *T 



(3) 



To, So and To are the cohesive parameters. They are all 

triaxiality dependent. Gy is the material yielding 
stress. Dcdi is an important parameter in a GTN 
simulation. It means the element height. The 
modification concerns the second part of equation (3), 
i.e. the calculation of So. A factor 1/5 is introduced 
because of convergence problems. In reference 
(Scheider 2009) the examined material is S460N, the 

yield stress is Gy = 470MPa and the element height is 
Dcdi = 0.1mm. Also H has a range boundary (1 < H< 4). 
If the triaxiality value jumps out of the range, the 
cohesive parameters are to be calculated under 
conditions of H = 1 or H = 4, respectively. 

Numerical Model 

Overview 

The conventional TSLs cannot describe the cyclic 
degradation for the cohesive element. In a trial to 
overcome this defect, in this paper, a tentative idea 
from Roe and Siegmund (2003) is taken as a starting 
point. Consulting the damage theory of the continuum 
element, a damage variable and a corresponding 
damage evolution law are introduced in the CZM 
following ideas of Scheider et al. (2003). This cyclic 
CZM is applied for fatigue analysis. 

There are two main fatigue phenomena - fatigue crack 
initiation and fatigue crack growth. In engineering 



practice, fatigue crack initiation is usually investigated 
using smooth unnotched specimens and fatigue crack 
growth is measured using pre-cracked specimens. The 
triaxiality conditions for these two standard specimens 
are totally different. Realistic modelling of the fatigue 
phenomena therefore requires taking the triaxiality 
into account in the damage evolution process. 

A triaxiality dependent cyclic CZM is proposed in this 
paper. The model should be used to reproduce the 
experimental fatigue phenomena which contain 
fatigue crack initiation and fatigue crack growth, just 
using a unique set of parameters. For simplicity, in the 
proposed CZM only mode-I fracture is considered. 
This means in CZM, just normal separation occurs. 

Construction of the Cyclic CZM 

With the damage variable D, the maximum cohesive 
traction and normal stiffness of cohesive elements will 
degrade. 



1 



h (H) 



(1-D) 



k N= k N (l-D) = ^(l-D) = 2 ^- 



T 



(4) 



(5) 



and *° are the initial and current maximum 
cohesive traction. The subscript (H) represents that the 
maximum cohesive traction is dependent on the 

k k 

triaxiality conditions. N and N are the initial and 

current normal stiffness of the cohesive element. 1 is 
a shape parameter for the TSL and relates to the 
critical separation So. 

For the triaxiality dependent cyclic CZM, it is difficult 
to judge the loading and unloading conditions because 
the separation S and the critical separation So alter 

S N 

simultaneously. So, a relative normal separation re 

AS N 



and its increment 



S N = 



are defined as 



°0 (H) 



AS N = 



oN 
(H ( ) 



_ °re (0 °re (t-Af) 



(6) 



(7) 



is the absolute normal separation. is the critical 
normal separation. The subscript (H) represents that 
the critical separation is triaxiality dependent. The 
subscripts t-At and t represent the previous and 
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S N S N 

current time increment. '~ M and ' are the absolute 
normal separation in the time increment t-At and t. 

1 and 1 '' are the critical normal separations 
in the time increment t-At and t. In the different time 
increments if the triaxiality value does not change, 
these two critical separations should be the same. 

<5 N S N 
re v' M ) and ,e ' f ' are the relative normal separation in 

the time increment t-At and t. 

AS N 

Based on the definition of rc , the loading and 
unloading conditions are summarized as follow. A 



loading step is indicated by 

i j- ■ ■ j- , ju AJ N <0 

unloading step is indicated by " 



AS N >0 



and an 



The response of the cohesive element in the first 
loading step should comply with the monotonic TSL, 
which is defined in equation (1). The following 
unloading and reloading process use the same 
equation, a linear function 



(8) 



Z 



N 



'~ M and ' are the normal traction in the different 
time increment. 

In the previous loading process a maximum relative 

5* 



normal separation rem * x may be reached. During the 
reloading period, if the current relative normal 
separation exceeds this maximum value, the reloading 
equation uses a new definition 



V* 5 ! J 



V 5 l J 



S N >& 



(9) 

■* and ° are the normal traction and separation in 
the current calculation increment. 

During the unloading/compression process, the 
normal separation maybe becomes negative value. At 
this moment a penalized stiffness should be adopted. 
But in the current cohesive model, the initial stiffness 
in normal direction is very high. So, application of the 
penalized stiffness is not necessary. 

In order to determine the current damage state, a 
damage evolution law should be included in the cyclic 
CZM. Here, a new damage evolution equation is 
proposed as 



D = 



A H e 



H 



I T o N (H) ](l-D) I?V 



D>0 



(10) 



D is the damage evolution. It is always a positive 

value. e is the material endurance limit. H is the 
triaxiality value. A and m are the material dependent 
damage controlling parameters. A is the material 
dependent triaxiality influence coefficient. 

After the cohesive element is broken, in the unloading 
or compression period in order to avoid the adjacent 
continuum element penetrating each other, a contact 

condition should be considered. Here, when & < , 
the contact is defined as 



(11) 



The application of the triaxiality dependent behaviour 
in cyclic CZM is a crucial aspect. The triaxiality 
dependent behaviour obtained from GTN simulation 
is just proper for monotonic loading condition. If this 
behaviour is applied for cyclic loading , the first 
problem is how to consider the triaxiality influence in 
unloading and reloading process. Actually, in the 
unloading/reloading process, it is still an open topic 
how the triaxiality influences the cohesive parameters. 
So, in this paper, a relative simple method is adopted. 
The triaxiality value keeps constant in the unloading 
process and a part of the reloading process. This part 
of the reloading process is also defined by equation (9), 
when the current relative normal separation does not 

S N 

exceed the maximum value " max , which is obtained 
from the previous loading step. In this case the 

S N 

triaxiality does not change. Once " max is exceeded, 
the triaxiality updates again. 

A schematic illustration for the cyclic process is shown 
in Fig. 2. 




one cycle 



FIG. 2 SCHEMATIC ILLUSTRATION FOR CYCLIC PROCESS OF 
CZM 
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The new cyclic CZM is contained in a user defined 
subroutine, and it is connected with the commercial 
finite element program ABAQUS. The fatigue damage 
can be calculated for every time increment. 

Model Verification 

Material Behaviour and Damage Parameters Analysis 

The triaxiality dependent cyclic CZM has been 
proposed and this new model should be verifyied by 
numerical simulation. The simulations contain two 
parts - fatigue crack initiation and fatigue crack 
growth. The primary idea is to reproduce the 
experimental results in fatigue test. The material 
S460N is chosen for the study. The triaxiality 
dependent behaviour for the simulated material is 
taken from reference (Scheider 2009). Some basic 
mechanical properties of this material are shown in 
Table 1. 

Table 1 Mechanical Properties For Material S460n (Boller 
And Seeger 1987) 



1000 



Materi 
al 


E(MPa) 


(jy(MP 
a) 


R»(MP 
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K'(MP 
a) 


«' 


MP 

a) 


S460N 


208000 


470 


682 


1181 


0.161 


312 



E is the material's Young's modulus. Y is the material 
yield stress. R m is the material ultimate tensile strength. 
K' is the material cyclic hardening coefficient, n is the 
material cyclic hardening exponent. o> is the material 
endurance limit. 

Before applying the cyclic CZM, it is necessary to 
discuss the material behaviour of the continuum 
element. The finite element model consists of two 
plane strain elements connected by one cohesive 
element, called homogeneous model. The model is 
loaded by cyclic displacement, the loading ratio is jR— 
1. 

In order to reflect the real material response under 
cyclic loading, the material behavior for the 
continuum elements uses a nonlinear kinematic 
hardening model in ABAQUS. Three different ways 
are offered by ABAQUS to define the nonlinear 
kinematic hardening component, specifying the 
material parameters directly, specifying half-cycle test 
data or specifying test data from a stabilized cycle. In 
this paper, nonlinear kinematic hardening model is 
identified from the stabilized cyclic test. 

The stabilized hysteresis loops for the response of the 
continuum element under three loading levels are 
depicted in Fig.3. 
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FIG. 3 STABILIZED HYSTERESIS LOOPS OF ONE CONTINUUM 
ELEMENT 

For the chosen material, the unknown damage 
parameters in the damage evolution law are just A, m 
and A. These parameters should be determined from 
the corresponding experiment test. However, it is 
impossible to discuss the influence of three damage 
parameters together. Hence, the material dependent 
triaxiality influence coefficient A is fixed first to A=1.2. 
The determination of this parameter is based on many 
trial calculations. The suggested A value for material 
S460N ranges between 0.7 and 1.2, a change of the 
parameter A will influence the other two parameters A 
and m. After A is chosen, a simulation is applied to 
analyze the material dependent damage controlling 
parameters A and m. The calculation model is still the 
homogeneous model. 
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FIG. 4 EFFECT OF DAMAGE PARAMETERS A AND m IN 
DAMAGE EVOLUTION LAW 

Series of the parameters A and m are investigated and 
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the calculation results are plotted in the material 
strain-life curves, shown in Fig. 4. 

It is obvious that the parameter m influences the slope 
of the strain-life curve. With increasing the value of m, 
the slope of the strain-life curve becomes flatter. The 
parameter A does not influence the slope of the strain- 
life curve. However, it makes the curve moving 
parallelly. According to these rules, the parameters A 
and m can be identified by fitting the numerical results 
to the experimental strain-life curve. 

Fatigue Crack Initiation Simulation 

In practice, for fatigue crack initiation investigation 
usually an uncracked smooth specimen is used in the 
experiment. The detailed dimension and geometry of 
the specimen is depicted in Fig. 5. The surface of the 
specimen is mechanically polished. The loading 
condition in experiment employs total strain control 
and the strain ratio is -1. The failure criterion in 
experiment is crack initiation. When a 0.5mm length 
crack is detected in the surface of the middle section, 
the experiment is stopped. 

t'rtt tfttgth hftwifn 
the cksnps 

u 



fa a 



v 
I 



r> — 0- — 

_L M 



FIG. 5 SPECIMEN IN FATIGUE EXPERIMENT 

Table 2 Experimental Strain-Life Curve Parameters (Boller 
And Seeger 1987) 



Material 


E(MPa) 


ff/(MPa) 


£/' 


b 


c 


S460N 


208000 


1218 


0.452 


-0.104 


-0.536 



All the test results are plotted in Fig. 6 together with 
the regression curve 



fatigue ductility coefficient, c is the fatigue ductility 
exponent. Values are taken from reference (Boiler and 
Seeger 1987) and shown in Table 2. 

10' 




FIG. 6 EXPERIMENTAL STRAIN-LIFE CURVE FOR S460N 
(Vormwald 2011) 

The simulation model is used here as described in the 
previous section. Two plane strain elements are 
connected by one cohesive element. When the 
cohesive element is broken, it means the crack 
initiation. Now the investigations just focus on the 
very low cycle fatigue regime, so the calculations are 
applied within 600 cycles. In the previous section, the 
analysis of the damage parameters has been explained. 
First, A is fixed to A=1.2 By fitting the simulation 
results to the experimental strain-life curve, the 
damage controlling parameters A and m can be 
determined. For material S460N, A=25 and m=3.7. The 
experiment and simulation results are plotted in Fig. 7. 
These parameters need further validation. 



S460N 



regression from eq.(12) 
experimental data 




10 100 
number of cycles 

FIG. 7 COMPARISON OF STRAIN-LIFE CURVE FOR S460N 



e = e ,+s ,=^(2N) b +£'(2N) C 
" °' el E y ' fK ' (12) 

" is the strain amplitude. aM and a,pl are the elastic 
and plastic part for the strain amplitude. N is the 

fatigue crack initiation life. 1 is the fatigue strength 

coefficient, b is the fatigue strength exponent, 1 is the 



Fatigue Crack Propagation Simulation 

As an important part of fatigue fracture, the fatigue 
crack growth simulation is applied to validate the 
damage parameters. 

Fatigue crack growth is usually measured by pre- 
cracked specimens, such as CT specimen, MT 
specimen or DCB specimen. The triaxiality condition is 
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totally different from the fatigue crack initiation 
simulation. The triaxiality dependent cyclic CZM 
should be able to reflect this different mechanism 
successfully. 

The experimental data of fatigue crack propagation on 
material S460N are shown in Fig. 8. It is clear that the 
experimental data can be expressed by Paris law. Now 
all the investigations focus on very low cycle fatigue 
regime, so just the experimental data at fast crack 
growth rates are used for comparison. 



l.E-07 




1 

Ay e „(N/mm) 



FIG. 8 EXPERIMENTAL FATIGUE CRACK GROWTH RATES 
CURVE FOR S460N (VORMWALD 2011) 

The simulation model is a two dimensional compact 
tension specimen. The width of the model is 50mm 
and the height is 60mm. The initial crack length ao is 
25mm and the ligament length is 25mm, too. The type 
of the continuum element is plane strain element. The 
cohesive elements locate along the ligament and the 
size of the single element is 0.125mm. The details for 
the finite element model are depicted in Fig. 9. 



boundary condition 
u,=0 



crack 



boundary condition 
u,=u v =0 




cohesive elements 



FIG. 9 FINITE ELEMENT MODEL FOR SIMULATED COMPACT 
TENSION SPECIMEN 

In the simulation, three different cyclic loading ranges 
are used: AF=1200N, AF=950N and AF=800N. The 
loading ratio is R=0. 

The response of the cohesive element needs an 
explanation. The CT model (in Fig. 9) under loading 
range AF=1200N is discussed. The first cohesive 
element at the crack tip is selected. The traction 



separation responses in two integration points of this 
cohesive element are plotted in Fig 10 and Figll. In 
the following figures, the two integration points of the 
cohesive element are called IP1 and IP2. The point IP2 
is closer to the crack tip. 
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FIG. 10 TRACTION SEPARATION RESPONSE IN IP2 

2000 i 1 




0.005 0.01 0.015 0.02 0.025 
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FIG. 11 TRACTION SEPARATION RESPONSE INIP1 

The traction separation responses in two integration 
points are very similar. After maximum traction To in 
the cohesive element is reached, the cohesive traction 
in the later loading step always degrades. Until the 
final failure, the traction becomes zero. The stiffness 
degradation with the damage variable is not obvious 
because the initial stiffness of the cohesive element is 
very high. Until the final failure, the separation in IP2 
is much bigger than the separation in IP1. 

The damage evolution in the two integration points 
are showed in Fig. 12. 
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FIG. 12 DAMAGE EVOLUTION IN TWO INTEGRATION POINTS 
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The damage evolving in IP2 is faster than in IP1. This 
is correct because IP2 is closer to the crack tip. During 
the crack growth, the failure sequence is from IP2 to 
IP1. So, the damage level in IP2 is always higher. The 
other feature in Fig. 12 should also be paid attention. 
The damage accumulation in cohesive element has a 
ladder-shaped distribution. The reason is that the 
damage does not evolve in unloading and a part of 
reloading process. 

The fatigue crack growth rates curve is the key point 
for comparison. The fatigue crack growth rate Aa/AN 
can be calculated directly, but the cyclic A/-integral 
will be computed by another approach. A finite 
element model of CT specimen without cohesive zone 
is built. The crack length inserted in this model is 
(flo+Aa). For one cycle the cyclic energy integral AJ can 
be computed by equation (13) (Wiithrich 1982). 

dAu. 



A/= M AWdy-At. '-ds 



AW= f 
Jo 



Aa.dAe.. 



(13) 



Aa 



As.. 



7 and 9 are cyclic stress and strain range, 
respectively. AW is the cyclic strain energy density, x 
and y are the Cartesian coordinates with the x-axis 
parallel to the crack surface. F is the integration path. 

ds is the line element lying on the integration path. ^' 



is the cyclic stress vector on the integration path, 
is the cyclic displacement variation. 



Au, 
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FIG. 13 COMPARISON OF FATIGUE CRACK GROWTH RATES 
CURVE FOR S460N 

The simulation results under three loading ranges are 
plotted together with the scatter band of experimental 
data in Fig. 13. The experimental Paris law can be 
reproduced by the CZM simulation to some extent. 
The deviation between experiment and simulation can 
also be accepted. The simulated crack growth rates 
curve seems to enter into an unstable region of fatigue 



crack growth. 

The simulated fatigue crack growth is faster than the 
experiment. One possible reason for this deviation is 
from the influence of the model dimension. In the 
experiment, the specimen is a three dimensional 
model, but in the simulation a two dimensional model 
is used. The calculation of the A/-integral leads to 
larger values for two dimensional problems than for 
three dimensional problems. Besides, the CZM may 
express the real cyclic process inaccurately. In the 
future more effort will be put into the improvement of 
the model. 

Conclusion and Outlook 

As a tool for material monotonic fracture analysis, the 
CZM is developed for many years and now it becomes 
a mature model. But for fatigue investigation, the 
CZM is just in a beginning state. In this paper, taking 
an idea from the pioneer as reference, a cyclic CZM 
including a new damage evolution law is proposed. 
This CZM involves the triaxiality influence 
characteristic and can be applied for various triaxiality 
conditions. In the very low cycle fatigue regime, the 
triaxiality dependent cyclic CZM is employed for 
fatigue crack initiation and fatigue crack growth 
simulation, reasonable agreement between experiment 
and simulation can be obtained. 

The new cyclic CZM is a developing numerical model 
and many improvements can be done in the future. 
The format of the damage evolution law may be 
modified. The triaxiality considered presently has a 
range boundary 1<H<4. Using GTN model simulations 
this range can be extended. The triaxiality dependent 
behavior for unloading/ compression is still a blank 
research field. If using the method in this paper, the 
triaxiality value keeps constant in the unloading 
process and a part of the reloading process. The 
triaxiality in the whole cyclic process will be 
discontinuous. This discontinuous triaxiality will lead 
to convergence problems. The triaxiality dependence 
characteristic for cohesive parameters is just taken 
from GTN model simulation. In reality, the triaxiality 
influence for cohesive parameters may be weaker. So, 
the correlative experiments need to be carried out in 
the future. 
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